
!Long-Term Care Insurance and the Family
!Corina Mommaerts

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
program conssav
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

use params
use wage        
use grids
use predata
use simulate

implicit none

include 'nlopt.f'
include 'mpif.h'

integer*8 opt
integer ires
real(prec) nlopt_tol, minf, objval
real(prec) x(Nparest) , x2(Nparest)
real(prec) lb(Nparest), ub(Nparest)

!for identification graphs:
integer iden !counter
integer, parameter :: niden=5
real(prec) identification_x1(niden), identification_x2(niden), identification_x3(niden)
real(prec) identification_y1(niden), identification_y2(niden), identification_y3(niden) 
real(prec) identification_y4(niden), identification_y5(niden), identification_y6(niden)

!for standard errors:
integer xval
real(prec) diffstep(Nparest)
real(prec) moments_plus(Nmoments), moments_minus(Nmoments)
real(prec) gradient_matrix(Nmoments,Nparest+1), gradient_matrix_transpose(Nparest+1,Nmoments)
real(prec) outer_a(Nparest+1,Nmoments),outer_b(Nparest+1,Nparest+1),outer(Nparest+1,Nparest+1),middle(Nparest+1,Nparest+1)
real(prec) varcov_parests(Nparest+1,Nparest+1)
real(prec) lambda(Nparest+1,Nmoments)
integer ErrorFlag

integer NTHREADS, OMP_GET_NUM_THREADS,OMP_GET_THREAD_NUM
integer iprov

integer algorithm
integer iseed

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!allocate arrays
 call arrays !in params.f90
 call allocate_wageproc(ierr)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
if (parallel.eq.1) then 
  call MPI_INIT(ierr) 							!initialize parallel
  call MPI_COMM_RANK(MPI_COMM_WORLD,pid,ierr)   !assign processor IDs (pid)
  call MPI_COMM_SIZE(MPI_COMM_WORLD,nproc,ierr) !determine number of processors started

  !$OMP PARALLEL
  NTHREADS = OMP_GET_NUM_THREADS()
  PRINT *, 'node',pid,'thread',OMP_GET_THREAD_NUM(),'of',NTHREADS
  !$OMP END PARALLEL

  if (nproc.ne.9 .and. nproc.ne.3 .and. nproc.ne.2 .and. nproc.ne.5 .and. nproc.ne.4) then
    node1 = 1
    node2 = nbignodes
    simnode1=1
    simnode2=nosims
  else if (nproc.eq.2) then
    if (pid.eq.0) then 
      node1 = 1
      node2 = 5
      simnode1=1
      simnode2=5000
    else
      node1 = 6
      node2 = 9
      simnode1=5001
      simnode2=9999
    endif
  else if (nproc.eq.4) then
	if (pid.eq.0) then 
      node1 = 1
      node2 = 3
      simnode1=1
      simnode2=2499
    else
	  node1 = 2*pid+2
	  node2 = 2*pid+3
	  simnode1 = 2500*pid
	  simnode2 = 2500*(pid+1)-1
	endif	
  else if (nproc.eq.5) then
	if (pid.eq.0) then 
      node1 = 1
      node2 = 1
      simnode1=1
      simnode2=1999
    else
	  node1 = 2*pid
	  node2 = 2*pid+1
	  simnode1 = 2000*pid
	  simnode2 = 2000*(pid+1)-1
	endif
  else
    node1=(nbignodes/nproc)*pid + 1
    node2=(nbignodes/nproc)*(pid+1)
    simnode1=(nosims/nproc)*pid+1
    simnode2=(nosims/nproc)*(pid+1)
  endif

  print*,'Process',pid,"of",nproc,'node range:',node1,node2
else
  pid=0
  node1=1
  node2=nh
  nproc=1
  print*,'not parallelizing.'
endif
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!


!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!read in inputs (in predata.f90)
call readin
call readin_moments
call readin_varcov

!set random numbers for simulations
!1. income shocks
seed = 1000
call ran1(seed,randincktype,nosims)
seed = 1001
call ran1(seed,randincp,nosims*(nt))
seed = 1002
call ran1(seed,randinck,nosims*(nt))
!2. survival rates
seed = 1003
call ran1(seed,randsurv,nosims*(nt))
!3. health shocks
seed = 1004
call ran1(seed,randhlth,nosims*(nt))
!4. love shocks
seed = 1005
call ran1(seed,randlovep,nosims*(nt))
seed = 1006
call ran1(seed,randlovek,nosims*(nt))
!5. HRS person profile random shock
seed = 1007 
call ran1(seed,randprofile,nosims*(nt))
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

whichiter=1

!create life expectancy table (using simulated people) - in pre-data.f90
call life_expectancy_table

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!PARAMETERS:
                                                                     
! baseline estimates
x(1)= -2E-7 					!formal care preference                     
x(2)= -3E-7 					!stigma
x(3)= 7.86604573807630075E-002  !altruism
x(4)= 2.26 				        !alpha
x(5)= 0.54142202266996531       !L shifter
x(6)= 1.8888369315929308        !gammab
x(7)= 2.1143038768435756        !gamma
x(8)= 23381				        !cminp
x(9)= -4.45076144931274663E-006 !mcaidpref
musimstartparam= 7 				!switch to 9 for joint identification (pareto grid modified) 								                                               				   	  

! ! alternative: parsimonious but just setting cminp (end result)
! x(1)= -4.99999998737621354E-007  !formal care preference                     
! x(2)= -2.31755258884958142E-007  !stigma
! x(3)= 6.52909264608252521E-002   !altruism
! x(4)= 2.2454939875410576 		   !alpha
! x(5)= 0.52208026120796114        !L shifter
! x(6)= 1.9153242605260958         !gammab
! x(7)= 2.1499999999999999 		   !gamma
! x(8)= 7124 					   !cminp
! x(9)= 2.67532961399935051E-005   !mcaidpref
! musimstartparam=7   

! ! alternative: parsimonious but just setting curvature on leisure
! x(1)= -9.22702436298188783E-008  !formal care preference                     
! x(2)= -3.95700813066076525E-007  !stigma
! x(3)= 8.83297278521432694E-002   !altruism
! x(4)= 2.22 				       !alpha
! x(5)= 0.52452760194464398        !L shifter
! x(6)=1.8970336079047210          !gammab
! x(7)=2.1003614151896874          !gamma
! x(8)=23580.293461533460		   !cminp
! x(9)= -9.85419745811059997E-006  !mcaidpref
! musimstartparam=7 
                                          

! ! alternative: parsimonious but just setting altruism
! x(1)= -3E-7  						!formal care preference                     
! x(2)= -3E-7  						!stigma
! x(3)= 0.04  						!altruism
! x(4)= 2.2716319329193899     		!alpha
! x(5)= 0.51      					!L shifter
! x(6)=1.905    					!gammab
! x(7)=2.1143038768435756         	!gamma
! x(8)=23381				        !cminp
! x(9)=-5E-6 						!mcaidpref
! musimstartparam=7 

! ! alternative: parsimonious but just setting pareto weight
! x(1)= -1.50649254645401872E-007 !formal care preference                     
! x(2)= -2.35491161679133947E-007 !stigma
! x(3)= 0.10524334334201949       !altruism
! x(4)= 2.2649312196221461        !alpha
! x(5)= 0.56242948243875890       !L shifter
! x(6)= 1.8732467274605571        !gammab
! x(7)= 2.1426693815587670        !gamma
! x(8)= 22465.366837295674        !cminp
! x(9)=-5.96653245991259973E-006  !mcaidpref
! musimstartparam=6
                 
! ! alternative: parsimonious but just setting gamma
! x(1)= -3.30352214333573249E-007 !formal care preference                     
! x(2)= -1.14966212773440131E-007 !stigma
! x(3)= 0.14999999999999999 	  !altruism
! x(4)= 2.2925542126723339  	  !alpha
! x(5)= 5.32578864270207911E-002  !L shifter
! x(6)= 1.7221220789081018        !gammab
! x(7)= 3.0  					  !gamma
! x(8)= 22277.357351495815 		  !cminp
! x(9)= 6.52323798844516695E-005  !mcaidpref
! musimstartparam=7 
				 
! ! alternative: parsimonious but NOT setting gamma
! x(1)= -4.99999998737621E-07 !formal care preference                     
! x(2)= -3.9198786470944E-07  !stigma
! x(3)= 0.04  				  !altruism
! x(4)= 2.22        		  !alpha
! x(5)= 0.344715581522315     !L shifter
! x(6)= 1.95        		  !gammab
! x(7)= 2.18899930607583      !gamma
! x(8)= 7124 !23381			  !cminp
! x(9)=0.0000166955971889051  !mcaidpref
! musimstartparam=6

! ! alternative: parsimonious (including gamma!)
! x(1)= -9.44201158602566285E-008 !formal care preference                     
! x(2)= -4.99999998737621354E-007 !stigma
! x(3)= 3.99999991059303284E-002  !altruism
! x(4)= 2.2200000286102295        !alpha
! x(5)= 0.46765124111122219       !L shifter
! x(6)= 1.7074910258068248        !gammab
! x(7)= 3.0         			  !gamma
! x(8)= 7124 			          !cminp
! x(9)= 4.79377556352439434E-005  !mcaidpref
! musimstartparam=6

if (simlowfloor.eq.1) then
	x(8)=0.1*x(8)
	print*,'LOW FLOOR:',x(8)
endif

if (annualperiods.eq.1) then
	x(1)=x(1)/2
	x(2)=x(2)/2
	x(8)=x(8)/2
	x(9)=x(9)/2
endif

diffstep=(/1E-8,1E-8,1E-3,.1,.01,0.1,0.1,1000.0,1E-7/)

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!


!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!create moment vector and weight matrix
weight_matrix=0.0
do momi=1,9
	moment_data_vector(   momi)=kasset_c1_p50_data(momi,1)
	moment_data_vector(9 +momi)=kasset_c2_p50_data(momi,1)
	moment_data_vector(18+momi)=kasset_c3_p50_data(momi,1)
	weight_matrix(   momi,   momi)=kasset_c1_p50_data(momi,2)
	weight_matrix(9 +momi,9 +momi)=kasset_c2_p50_data(momi,2)
	weight_matrix(18+momi,18+momi)=kasset_c3_p50_data(momi,2)
enddo
do momi=1,7
	moment_data_vector(27+momi)=kasset_c4_p50_data(momi,1)
	weight_matrix(27+momi,27+momi)=kasset_c4_p50_data(momi,2)	
enddo
do momi=1,5
	moment_data_vector(34+momi)=kasset_c5_p50_data(momi,1)
	weight_matrix(34+momi,34+momi)=kasset_c5_p50_data(momi,2)
enddo
do momi=1,9
	moment_data_vector(39+momi)=passet_c1_p50_data(momi,1)
	moment_data_vector(48+momi)=passet_c2_p50_data(momi,1)
	moment_data_vector(57+momi)=passet_c3_p50_data(momi,1)
	weight_matrix(39+momi,39+momi)=passet_c1_p50_data(momi,2)
	weight_matrix(48+momi,48+momi)=passet_c2_p50_data(momi,2)
	weight_matrix(57+momi,57+momi)=passet_c3_p50_data(momi,2)
enddo
do momi=1,7
	moment_data_vector(66+momi)=passet_c4_p50_data(momi,1)
	weight_matrix(66+momi,66+momi)=passet_c4_p50_data(momi,2)	
enddo   
do momi=1,5
	moment_data_vector(73+momi)=passet_c5_p50_data(momi,1)
	weight_matrix(73+momi,73+momi)=passet_c5_p50_data(momi,2)	
enddo 
!these next numbers come from stata do files
moment_data_vector(79)=.0269316  !insurance -1
moment_data_vector(80)=.040744   !insurance -2
moment_data_vector(81)=.0598065  !insurance -3
moment_data_vector(82)=.0863656  !insurance -4
moment_data_vector(83)=.1699518  !insurance -5
moment_data_vector(84)=.3927813  !informal care - 1
moment_data_vector(85)=.4819159  !informal care - 2
moment_data_vector(86)=.6164384  !informal care - 3
moment_data_vector(87)=.5985915  !informal care - 4
moment_data_vector(88)=.5391621  !informal care - 5
moment_data_vector(89)=.2031141  !Medicaid
moment_data_vector(90)=.6874669  !Kid FT - cohort 1
moment_data_vector(91)=.6714159  !Kid FT - cohort 2
moment_data_vector(92)=.5942078  !Kid FT - cohort 3
moment_data_vector(93)=.5285962  !Kid FT - cohort 4
moment_data_vector(94)=.4501608  !Kid FT - cohort 5

weight_matrix(79,79)=33081.271  !insurance -1
weight_matrix(80,80)=32979.033  !insurance -2
weight_matrix(81,81)=33212.72   !insurance -3
weight_matrix(82,82)=33314.958  !insurance -4
weight_matrix(83,83)=33344.168  !insurance -5
weight_matrix(84,84)=5824.0171  !informal care -1
weight_matrix(85,85)=4216.5389  !informal care -2
weight_matrix(86,86)=3008.8694  !informal care -3
weight_matrix(87,87)=2341.1477  !informal care -4
weight_matrix(88,88)=2262.8347  !informal care -5
weight_matrix(89,89)=71018.7    !Medicaid
weight_matrix(90,90)=8081.0668  !Kid FT - cohort 1
weight_matrix(91,91)=9598.1365  !Kid FT - cohort 2
weight_matrix(92,92)=12542.534  !Kid FT - cohort 3
weight_matrix(93,93)=9863.0895  !Kid FT - cohort 4
weight_matrix(94,94)=5316.154   !Kid FT - cohort 5

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!


!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
if (DoEstimation.eq.1) then

  iseed=7121986
  call nlosr(iseed)

  call nlo_create(opt, NLOPT_GN_CRS2_LM, Nparest)  !global optimizer
  call nlo_get_algorithm(algorithm, opt)

  call nlo_get_lower_bounds(ires,opt,lb)
  call nlo_get_upper_bounds(ires,opt,ub)
  lb(1)= -5E-7    !formal care preference
  ub(1)= 0.0d0    !formal care preference
  lb(2)= -5E-7    !stigma
  ub(2)= 0.0d0    !stigma
  lb(3)= 0.01d0   !altruism
  ub(3)= 0.15d0   !altruism
  lb(4)= 2.1d0 	  !alpha
  ub(4)= 2.3      !alpha
  lb(5)= 0.45d0   !L shifter
  ub(5)= 0.65d0   !L shifter
  lb(6)= 1.6d0 	  !gammab
  ub(6)= 2.0d0    !gammab
  lb(7)= 1.95d0   !gamma
  ub(7)= 2.25d0   !gamma
  lb(8)= 15000d0  !cminp
  ub(8)= 25000d0  !cminp
  lb(9)= -7E-5    !mcaidpref
  ub(9)= -1E-7    !mcaidpref 
  if (estimatesimple==1) then 
	ub(1)= 8E-7   !formal care preference
	lb(5)= 0.05d0 !L shifter
	lb(7)= 2.05d0 !gamma
	ub(7)= 3.05d0 !gamma
	lb(8)= 0d0 	  !cminp
	ub(9)= 1E-4   !mcaidpref 
  endif
  call nlo_set_lower_bounds(ires,opt,lb)
  call nlo_set_upper_bounds(ires,opt,ub)
  call nlo_set_min_objective(ires, opt, runmyfunction,1) !specify objective function

  nlopt_tol= 1.D-8 !for global optimizer
  call nlo_set_stopval(ires, opt, nlopt_tol)  !specify a stopping criteria
  call nlo_set_xtol_rel(ires, opt, nlopt_tol) !specify a stopping criteria
  call nlo_set_xtol_abs(ires, opt, nlopt_tol) !specify a stopping criteria
  call nlo_set_ftol_rel(ires, opt, nlopt_tol) !specify a stopping criteria
  call nlo_set_ftol_abs(ires, opt, nlopt_tol) !specify a stopping criteria

  !optimize! x is initial guess vector, then returns optimized values
  call nlo_optimize(ires,opt,x,minf)!ires>0 if success
  if (ires.lt.0) then
    write(*,*) 'nlopt failed!'
  else
    write(*,*) 'found min at ', x
    write(*,*) 'min val = ', minf
    write(*,*) 'return:',ires
  endif

else 
  call runmyfunction(objval,Nparest,x,0d0,0)
  
  if (DoIdentification==1) then
      
	!for mcaid preference
    identification_y1(1)=avinf(1)
    identification_y2(1)=avinf(2)	
    identification_y3(1)=avinf(3)
    identification_y4(1)=avinf(4)
    identification_y5(1)=avinf(5)
    identification_x1=(/x(9), x(9)-3D-006, x(9)-6D-006, x(9)+3D-006, x(9)+6D-006/)
	do iden=2,niden
		x(9)=identification_x1(iden)
		call runmyfunction(objval,Nparest,x,0d0,0)
		identification_y1(iden)=avinf(1)
		identification_y2(iden)=avinf(2)
		identification_y3(iden)=avinf(3)
		identification_y4(iden)=avinf(4)      
		identification_y5(iden)=avinf(5)      
	enddo
	if (pid==0) then 
		open(1,file='output/identification_mcaidpref.csv',action="write",status="replace")
		write(1,"(6(f0.12,',',:))") identification_x1(:)
		write(1,"(6(f0.12,',',:))") identification_y1(:)
		write(1,"(6(f0.12,',',:))") identification_y2(:)
		write(1,"(6(f0.12,',',:))") identification_y3(:)
		write(1,"(6(f0.12,',',:))") identification_y4(:)  
		write(1,"(6(f0.12,',',:))") identification_y5(:)  
	endif 
	
	! !for formal care preference
    ! identification_y1(1)=avinf(1)
    ! identification_y2(1)=avinf(2)	
    ! identification_y3(1)=avinf(3)
    ! identification_y4(1)=avinf(4)
    ! identification_y5(1)=avinf(5)
    ! identification_y6(1)=avinf(6)
    ! identification_x1=(/x(1), x(1)-2D-007, x(1)-3D-007, x(1)+2D-007, x(1)+4D-007/)
	! do iden=2,niden
		! x(1)=identification_x1(iden)
		! call runmyfunction(objval,Nparest,x,0d0,0)
		! identification_y1(iden)=avinf(1)
		! identification_y2(iden)=avinf(2)
		! identification_y3(iden)=avinf(3)
		! identification_y4(iden)=avinf(4)      
		! identification_y5(iden)=avinf(5)      
		! identification_y6(iden)=avinf(6)      
	! enddo
	! if (pid==0) then 
		! open(1,file='output/identification_formalpref.csv',action="write",status="replace")
		! write(1,"(6(f0.12,',',:))") identification_x1(:)
		! write(1,"(6(f0.12,',',:))") identification_y1(:)
		! write(1,"(6(f0.12,',',:))") identification_y2(:)
		! write(1,"(6(f0.12,',',:))") identification_y3(:)
		! write(1,"(6(f0.12,',',:))") identification_y4(:)  
		! write(1,"(6(f0.12,',',:))") identification_y5(:)  
		! write(1,"(6(f0.12,',',:))") identification_y6(:)  
	! endif 	  

	! !for stigma
    ! identification_y1(1)=avins(1)
    ! identification_y2(1)=avins(2)	
    ! identification_y3(1)=avins(3)
    ! identification_y4(1)=avins(4)
    ! identification_y5(1)=avins(5)
    ! identification_y6(1)=avins(6)
    ! identification_x1=(/x(2), x(2)-4D-007, x(2)-6D-007,  x(2)+3D-007 , x(2)+6D-007/)
	! do iden=2,niden
		! x(2)=identification_x1(iden)
		! call runmyfunction(objval,Nparest,x,0d0,0)
		! identification_y1(iden)=avins(1)
		! identification_y2(iden)=avins(2)
		! identification_y3(iden)=avins(3)
		! identification_y4(iden)=avins(4)      
		! identification_y5(iden)=avins(5)      
		! identification_y6(iden)=avins(6)      
	! enddo
	! if (pid==0) then 
		! open(1,file='output/identification_stigma.csv',action="write",status="replace")
		! write(1,"(6(f0.12,',',:))") identification_x1(:)
		! write(1,"(6(f0.12,',',:))") identification_y1(:)
		! write(1,"(6(f0.12,',',:))") identification_y2(:)
		! write(1,"(6(f0.12,',',:))") identification_y3(:)
		! write(1,"(6(f0.12,',',:))") identification_y4(:)  
		! write(1,"(6(f0.12,',',:))") identification_y5(:)  
		! write(1,"(6(f0.12,',',:))") identification_y6(:)  
	! endif 	

	! !for altruism
    ! ! identification_y1(1)=wealthcp50(3,3)
    ! ! identification_y2(1)=wealthck50(3,3)	
    ! ! identification_y3(1)=wealthcp50(1,7)
    ! ! identification_y4(1)=wealthck50(1,7)
    ! ! identification_y5(1)=wealthcp50(2,5)
    ! ! identification_y6(1)=wealthck50(2,5)
    ! ! identification_x1=(/x(3), 0.0D0, x(3)-2D-002, x(3)+1D-002, 0.15D0/)
	! ! do iden=2,niden
		! ! x(3)=identification_x1(iden)
		! ! call runmyfunction(objval,Nparest,x,0d0,0)
		! ! identification_y1(iden)=wealthcp50(3,3)
		! ! identification_y2(iden)=wealthck50(3,3)
		! ! identification_y3(iden)=wealthcp50(1,7)
		! ! identification_y4(iden)=wealthck50(1,7)     
		! ! identification_y5(iden)=wealthcp50(2,5)    
		! ! identification_y6(iden)=wealthck50(2,5)      
	! ! enddo
	! identification_y1(1)=wealthcp50(1,1)/wealthck50(1,1)
    ! identification_y2(1)=wealthcp50(1,2)/wealthck50(1,2)
    ! identification_y3(1)=wealthcp50(1,3)/wealthck50(1,3)
    ! identification_y4(1)=wealthcp50(2,1)/wealthck50(2,1)
    ! identification_y5(1)=wealthcp50(2,2)/wealthck50(2,2)
    ! identification_x1=(/x(3), 0.0D0, x(3)-2D-002, x(3)+1D-002, 0.15D0/)
	! do iden=2,niden
		! x(3)=identification_x1(iden)
		! call runmyfunction(objval,Nparest,x,0d0,0)
		! identification_y1(iden)=wealthcp50(1,1)/wealthck50(1,1)
		! identification_y2(iden)=wealthcp50(1,2)/wealthck50(1,2)
		! identification_y3(iden)=wealthcp50(1,3)/wealthck50(1,3)
		! identification_y4(iden)=wealthcp50(2,1)/wealthck50(2,1)
		! identification_y5(iden)=wealthcp50(2,2)/wealthck50(2,2)   
	! enddo
	! if (pid==0) then 
		! !open(1,file='output/identification_altruism.csv',action="write",status="replace")
		! open(1,file='output/identification_altruism_ratiomoments.csv',action="write",status="replace")
		! write(1,"(6(f0.12,',',:))") identification_x1(:)
		! write(1,"(6(f0.12,',',:))") identification_y1(:)
		! write(1,"(6(f0.12,',',:))") identification_y2(:)
		! write(1,"(6(f0.12,',',:))") identification_y3(:)
		! write(1,"(6(f0.12,',',:))") identification_y4(:)  
		! write(1,"(6(f0.12,',',:))") identification_y5(:)  
		! write(1,"(6(f0.12,',',:))") identification_y6(:)  
	! endif 
	
	! !for initial pareto weight
    ! identification_y1(1)=wealthcp50(1,1)/wealthck50(1,1)
    ! identification_y2(1)=wealthcp50(1,2)/wealthck50(1,2)
    ! identification_y3(1)=wealthcp50(1,3)/wealthck50(1,3)
    ! identification_y4(1)=wealthcp50(2,1)/wealthck50(2,1)
    ! identification_y5(1)=wealthcp50(2,2)/wealthck50(2,2)
    ! identification_x1=(/7,5,6,8,9/)
	! do iden=2,niden
		! musimstartparam=identification_x1(iden)
		! call runmyfunction(objval,Nparest,x,0d0,0)
		! identification_y1(iden)=wealthcp50(1,1)/wealthck50(1,1)
		! identification_y2(iden)=wealthcp50(1,2)/wealthck50(1,2)
		! identification_y3(iden)=wealthcp50(1,3)/wealthck50(1,3)
		! identification_y4(iden)=wealthcp50(2,1)/wealthck50(2,1)
		! identification_y5(iden)=wealthcp50(2,2)/wealthck50(2,2)   
	! enddo
	! ! identification_y1(1)=wealthcp50(3,3)
    ! ! identification_y2(1)=wealthck50(3,3)	
    ! ! identification_y3(1)=wealthcp50(1,7)
    ! ! identification_y4(1)=wealthck50(1,7)
    ! ! identification_y5(1)=wealthcp50(2,5)
    ! ! identification_y6(1)=wealthck50(2,5)
    ! ! identification_x1=(/7,5,6,8,9/)
	! ! do iden=2,niden
		! ! x(3)=identification_x1(iden)
		! ! call runmyfunction(objval,Nparest,x,0d0,0)
		! ! identification_y1(iden)=wealthcp50(3,3)
		! ! identification_y2(iden)=wealthck50(3,3)
		! ! identification_y3(iden)=wealthcp50(1,7)
		! ! identification_y4(iden)=wealthck50(1,7)     
		! ! identification_y5(iden)=wealthcp50(2,5)    
		! ! identification_y6(iden)=wealthck50(2,5)      
	! ! enddo
	! if (pid==0) then 
		! open(1,file='output/identification_pareto.csv',action="write",status="replace")
		! !open(1,file='output/identification_pareto_wealthmoments.csv',action="write",status="replace")
		! write(1,"(6(f0.12,',',:))") identification_x1(:)
		! write(1,"(6(f0.12,',',:))") identification_y1(:)
		! write(1,"(6(f0.12,',',:))") identification_y2(:)
		! write(1,"(6(f0.12,',',:))") identification_y3(:)
		! write(1,"(6(f0.12,',',:))") identification_y4(:)  
		! write(1,"(6(f0.12,',',:))") identification_y5(:)  
	! endif 

	! !for initial pareto weight+altruism constant part
    ! identification_y1(1)=wealthcp50(1,2)/wealthck50(1,2)
    ! identification_y2(1)=wealthcp50(1,3)/wealthck50(1,3)
    ! identification_y3(1)=wealthcp50(2,2)/wealthck50(2,2)
	! identification_y4(1)=wealthcp50(3,3)
	! identification_y5(1)=wealthcp50(1,7)
    ! identification_y6(1)=wealthcp50(2,5)
    ! identification_x1=(/9,7,8,10,11/)
	! do iden=2,niden
		! musimstartparam=identification_x1(iden) !new pareto weight
		! x(3)= (paretoweight(musimstartparam)*(1+0.6*0.0786604573807630075)-0.6)/(paretoweight(musimstartparam)*0.6) !altruism to keep constant 
		! if (x(3)<0) x(3)=0.0d0
		! print*,musimstartparam,paretoweight(musimstartparam),x(3)
		! call runmyfunction(objval,Nparest,x,0d0,0)
		! identification_y1(iden)=wealthcp50(1,2)/wealthck50(1,2)
		! identification_y2(iden)=wealthcp50(1,3)/wealthck50(1,3)
		! identification_y3(iden)=wealthcp50(2,2)/wealthck50(2,2)   
		! identification_y4(iden)=wealthcp50(3,3)
		! identification_y5(iden)=wealthcp50(1,7)
		! identification_y6(iden)=wealthcp50(2,5)
	! enddo
	! if (pid==0) then 
		! open(1,file='output/identification_paretoaltruismjoint.csv',action="write",status="replace")
		! write(1,"(6(f0.12,',',:))") identification_x1(:)
		! write(1,"(6(f0.12,',',:))") identification_y1(:)
		! write(1,"(6(f0.12,',',:))") identification_y2(:)
		! write(1,"(6(f0.12,',',:))") identification_y3(:)
		! write(1,"(6(f0.12,',',:))") identification_y4(:)  
		! write(1,"(6(f0.12,',',:))") identification_y5(:)  
		! write(1,"(6(f0.12,',',:))") identification_y6(:)  
	! endif 
	
	
    endif !end identification part

    !get standard errors
    if (DoStandardErrors==1) then

      !cycle through all parameters:
      do xval=1,Nparest
        
        x2=x

        !plus
        x2(xval) = x(xval) + diffstep(xval)
        call runmyfunction(objval,Nparest,x2,0d0,0)
        moments_plus=moment_simu_vector

        !minus
        x2(xval) = x(xval) - diffstep(xval)
        call runmyfunction(objval,Nparest,x2,0d0,0)
        moments_minus=moment_simu_vector

        gradient_matrix(:,xval) = (moments_plus-moments_minus)/(2*diffstep(xval))

      enddo
      
      !pareto weight
      x2=x
      !plus
      musimstartparam=10 !8
      call runmyfunction(objval,Nparest,x2,0d0,0)
      moments_plus=moment_simu_vector
      !minus
      musimstartparam=8 !6
      call runmyfunction(objval,Nparest,x2,0d0,0)
      moments_minus=moment_simu_vector      
      !gradient
      gradient_matrix(:,Nparest+1) = (moments_plus-moments_minus)/(2*0.1)

      !done cycling through all parameters
      gradient_matrix_transpose=transpose(gradient_matrix)

      !variance of estimates:
      ! outer part: outer=(G'WG)^{-1}
      outer_a=matmul(gradient_matrix_transpose,weight_matrix)
      outer_b=matmul(outer_a,gradient_matrix)
      call FindInv(outer_b, outer, Nparest+1, ErrorFlag)
      ! middle part: middle=G'Wvarcov_matrixWG
      outer_a=matmul(outer_a,varcov_matrix)
      outer_a=matmul(outer_a,weight_matrix)
      middle=matmul(outer_a,gradient_matrix)
      ! put together:
      varcov_parests=matmul(outer,middle)
      varcov_parests=matmul(varcov_parests,outer)
	  
	  !sensitivity matrix: lambda=-(G'WG)^{-1}G'W but forgot the -
	  lambda=matmul(outer,gradient_matrix_transpose)
	  lambda=matmul(lambda,weight_matrix)	  

	  print*,'var-cov matrix:'
      do xval=1,Nparest+1
        print*,varcov_parests(xval,:)
      enddo
	  
	  print*,'sensitivity matrix:'
	  do xval=1,Nparest+1
	    print*,lambda(xval,:)
	  enddo
	  
    endif !end standard errors

endif !end estimation or not
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

call MPI_FINALIZE(ierr)

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
end program conssav
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!


subroutine runmyfunction(valreturn,n,x,grad,need_gradient)
!this is the objective function
!val returns result of the objective function
!x is array of optimization parameters (size n)

!INCLUDED MODULES
use grids           !creates asset grid
use valuefunction   !solves value function (consumption) at given asset/income node
use params          !defines parameters
use simulate        !simulates the model after it has been solved
use simwrite        !writes the results of the simulations
use states          !defines state variables for easy transport between modules
use utilityfs
use wage            !creates income grid (in wageproc)

implicit none
include 'mpif.h'

!optimization variables
integer n,need_gradient
double precision, intent(out) :: valreturn
double precision x(n),grad

!DECLARE VARIABLES
integer ier

integer t !counter over time/age
integer h !counter over health nodes
integer m !counter over pareto weight
integer i !counter over insurance

integer outcsv1
integer ins_d

real(prec) temp1, temp2

integer t1,t2,clock_rate,clock_max

integer nodek,wknode,aknode
integer proc, OMP_GET_THREAD_NUM

integer count

real(prec) x_transfer(Nparest)
real(prec) valreturnmat(nm)

call system_clock ( t1, clock_rate, clock_max )

x_transfer=0.0
if (pid==0) x_transfer=x
call MPI_ALLREDUCE(x_transfer,x,Nparest,MPI_DOUBLE_PRECISION,MPI_SUM,MPI_COMM_WORLD,ierr)

formalpref  = x(1) 
stigma      = x(2) 
altruism    = x(3)
alpha       = x(4)
lshifter	= x(5)
termphiL=0.5 !old holdover
termphiH=0.5 !old holdover
gammab		= x(6)
gamma = x(7) 
cminp       = x(8)
cmink = cminp
mcaidpref   = x(9) 

if (estimatesimple==1) then  !toggle these on/off
	!altruism=0.04
	!alpha=2.22
	gamma=3.0
	!cminp=7124
	!cmink=7124
endif

if (pid==0) print*,formalpref,stigma,altruism,alpha,lshifter,gammab,gamma,cminp,mcaidpref,musimstartparam

!*****************************************************!
!***********SOLVE THE MODEL***************************!
!*****************************************************!

valuep_padj=-9999999.0
valuep_kadj=-9999999.0
valuek_padj=-9999999.0
valuek_kadj=-9999999.0

valuep=0
valuek=0

valuek_deadparent=0

valuekNT_deadparent=0
valuekPT_deadparent=0
valuekFT_deadparent=0

conexpkNT_deadparent=0
conexpkPT_deadparent=0
conexpkFT_deadparent=0

valuep_divorced=0
valuek_divorced=0

valuepNT_divorced=0
valuepPT_divorced=0
valuepFT_divorced=0
valuekNT_divorced=0
valuekPT_divorced=0
valuekFT_divorced=0

conexppNT_divorced=0
conexppPT_divorced=0
conexppFT_divorced=0
conexpkNT_divorced=0
conexpkPT_divorced=0
conexpkFT_divorced=0

valuepNT0_married=0
valuepPT0_married=0
valuepFT0_married=0
valuepNT100_married=0
valuepPT100_married=0
valuepFT100_married=0

valuekNT0_married=0
valuekPT0_married=0
valuekFT0_married=0
valuekNT100_married=0
valuekPT100_married=0
valuekFT100_married=0

valueNT0_married=0
valuePT0_married=0
valueFT0_married=0
valueNT100_married=0
valuePT100_married=0
valueFT100_married=0

conexppNT0_married=0
conexppPT0_married=0
conexppFT0_married=0
conexppNT100_married=0
conexppPT100_married=0
conexppFT100_married=0

conexpkNT0_married=0
conexpkPT0_married=0
conexpkFT0_married=0
conexpkNT100_married=0
conexpkPT100_married=0
conexpkFT100_married=0

kappaNT0_married=0
kappaPT0_married=0
kappaFT0_married=0
kappaNT100_married=0
kappaPT100_married=0
kappaFT100_married=0


!GENERATE INCOME GRID
 call wageproc
 
!**** (1) LOOP OVER TIME ***************************************************************************!
do t=nt,1,-1

  if (pid==0) print*,'AGE',t,pid

  age   = t   !STATE VARIABLE (in states.f90)
  ageto = t+1 !STATE VARIABLE (in states.f90)

  !define asset grid bounds
  call assetbounds
  call assetgrid 

  !*********************************************************************************************!
  !GET VALUE WHEN PARENT IS DEAD
  !$OMP PARALLEL DO PRIVATE(nodek)
  do nodek=1,nonodesk    
    call solvevalfun_deadparent(nodek)
  enddo !end loop over kid asset nodes
  !$OMP END PARALLEL DO
  !*********************************************************************************************!
  
  !*********************************************************************************************! 
  !GET VALUE OF "DIVORCE"
  !parent and kid values still linked through parent altruism (just bequest...)
  !must solve jointly - grid search
  do bignode=node1,node2
    health=ceiling(real(bignode)/real(nwp)) !(1=healthy,2=semi-unhealthy,3=very-unhealthy)
    wpnode=bignode-(health-1)*nwp           !parent income
    do apnode=1,na

      !$OMP PARALLEL DO PRIVATE(nodek)
      do nodek=1,nonodesk    
        call solvevalfun_divorced(nodek)

        if (noncoop==1) then
          do i=1,ni
            valuep(age,i,bignode,:,apnode,nodek)=valuep_divorced(age,i,bignode,apnode,nodek)
            valuek(age,i,bignode,:,apnode,nodek)=valuek_divorced(age,i,bignode,apnode,nodek)
            consump(age,i,bignode,:,apnode,nodek)=conexpp_divorced(age,i,bignode,apnode,nodek)
            consumk(age,i,bignode,:,apnode,nodek)=conexpk_divorced(age,i,bignode,apnode,nodek)
            LFP(age,i,bignode,:,apnode,nodek)=LFP_divorced(age,i,bignode,apnode,nodek)
            family(age,i,bignode,:,apnode,nodek)=0
            married(age,i,bignode,:,apnode,nodek)=20     
          enddo  
          if (age==1) insurance(bignode,:,apnode,nodek)=insurance_divorced(bignode,apnode,nodek)
        endif
      enddo !end loop over kid asset nodes
      !$OMP END PARALLEL DO

    enddo !end loop over parent asset nodes
  enddo !end loop over health / parent income / kid income states  
  !*********************************************************************************************!
  
  if (noncoop==1)  goto 1234

  !*********************************************************************************************!
  ! GET VALUE OF "MARRIAGE"
  do bignode=node1,node2

    health=ceiling(real(bignode)/real(nwp)) !(1=healthy,2=semi-unhealthy,3=very-unhealthy)
    wpnode=bignode-(health-1)*nwp         !parent income

    !**** LOOP OVER PARETO WEIGHT *****************************************************!
    !     (first time through to get value of marriage)
    do m=1,nm

      mnode=m
      mweight=paretoweight(m)

      do apnode=1,na

        !$OMP PARALLEL DO PRIVATE(nodek)
        do nodek=1,nonodesk
          call solvevalfun_married(nodek)
        enddo !end loop over kid asset
        !$OMP END PARALLEL DO
		
      enddo !end loop over parent asset
    enddo !end loop over pareto weight
    !**********************************************************************************!

    !**** LOOP OVER PARETO WEIGHT *****************************************************!
    !     (second time through for weight adjustments)
    do m=1,nm

      mnode=m
      mweight=paretoweight(m)

      do apnode=1,na

        do nodek=1,nonodesk
            wknode=ceiling(real(nodek)/real(na))
            aknode=nodek-(wknode-1)*na 

          !*** loop over insurance states ***!
          do i=1,ni

            if (age .eq. 1) then
              ins_d=insurance_divorced(bignode,apnode,nodek)
            else
              ins_d=i
            endif

            if (limcom.eq.0) then
              value(age,i,bignode,m,apnode,nodek)=value_married(age,i,bignode,m,apnode,nodek)
              valuep(age,i,bignode,m,apnode,nodek)=valuep_married(age,i,bignode,m,apnode,nodek)
              valuek(age,i,bignode,m,apnode,nodek)=valuek_married(age,i,bignode,m,apnode,nodek)
              consump(age,i,bignode,m,apnode,nodek)=conexpp_married(age,i,bignode,m,apnode,nodek)  
              consumk(age,i,bignode,m,apnode,nodek)=conexpk_married(age,i,bignode,m,apnode,nodek)  
              kappa(age,i,bignode,m,apnode,nodek)=kappa_married(age,i,bignode,m,apnode,nodek)      
              mu(age,i,bignode,m,apnode,nodek)=mweight          
              LFP(age,i,bignode,m,apnode,nodek)=LFP_married(age,i,bignode,m,apnode,nodek)
              family(age,i,bignode,m,apnode,nodek)=famcare_married(age,i,bignode,m,apnode,nodek)
              married(age,i,bignode,m,apnode,nodek)=10
            else
            !*** 3 SCENARIOS *************************!

            !!!!ONE: BOTH PREFER MARRIAGE --> STAY MARRIED
            if (valuep_married(age,i,bignode,m,apnode,nodek) .ge. valuep_divorced(age,ins_d,bignode,apnode,nodek) .and. & 
                valuek_married(age,i,bignode,m,apnode,nodek) .ge. valuek_divorced(age,ins_d,bignode,apnode,nodek)) then
                
              value(age,i,bignode,m,apnode,nodek)=value_married(age,i,bignode,m,apnode,nodek)
              valuep(age,i,bignode,m,apnode,nodek)=valuep_married(age,i,bignode,m,apnode,nodek)
              valuek(age,i,bignode,m,apnode,nodek)=valuek_married(age,i,bignode,m,apnode,nodek)
              consump(age,i,bignode,m,apnode,nodek)=conexpp_married(age,i,bignode,m,apnode,nodek)  
              consumk(age,i,bignode,m,apnode,nodek)=conexpk_married(age,i,bignode,m,apnode,nodek)  
              kappa(age,i,bignode,m,apnode,nodek)=kappa_married(age,i,bignode,m,apnode,nodek)      
              mu(age,i,bignode,m,apnode,nodek)=mweight          
              LFP(age,i,bignode,m,apnode,nodek)=LFP_married(age,i,bignode,m,apnode,nodek)
              family(age,i,bignode,m,apnode,nodek)=famcare_married(age,i,bignode,m,apnode,nodek)
              married(age,i,bignode,m,apnode,nodek)=10

            !!!!TWO: KID PREFERS DIVORCE --> WEIGHT MODIFICATION
            elseif (valuek_divorced(age,ins_d,bignode,apnode,nodek).ge.valuek_married(age,i,bignode,m,apnode,nodek)) then
                
              !first check to see if adjusted value not already calculated (skip if already calculated)
              if (valuek_kadj(age,i,bignode,apnode,nodek)<-100.0) then
                !find weight for which valuek_married=valuek_divorced
                call find_valkwgt(valuek_divorced(age,ins_d,bignode,apnode,nodek),i,nodek)
              endif
              valuep(age,i,bignode,m,apnode,nodek)=valuep_kadj(age,i,bignode,apnode,nodek)
              valuek(age,i,bignode,m,apnode,nodek)=valuek_kadj(age,i,bignode,apnode,nodek)
              consump(age,i,bignode,m,apnode,nodek)=conexpp_kadj(age,i,bignode,apnode,nodek)
              consumk(age,i,bignode,m,apnode,nodek)=conexpk_kadj(age,i,bignode,apnode,nodek)
              kappa(age,i,bignode,m,apnode,nodek)=kappa_kadj(age,i,bignode,apnode,nodek)
              mu(age,i,bignode,m,apnode,nodek)=mu_kadj(age,i,bignode,apnode,nodek)
              LFP(age,i,bignode,m,apnode,nodek)=LFP_kadj(age,i,bignode,apnode,nodek)
              family(age,i,bignode,m,apnode,nodek)=famcare_kadj(age,i,bignode,apnode,nodek)
              married(age,i,bignode,m,apnode,nodek)=11
              value(age,i,bignode,m,apnode,nodek)=mu(age,i,bignode,m,apnode,nodek)* &
                                                          valuep(age,i,bignode,m,apnode,nodek) + &
                                                        (1.0-mu(age,i,bignode,m,apnode,nodek))* &
                                                          valuek(age,i,bignode,m,apnode,nodek)

            !!!!THREE: PARENT PREFERS DIVORCE --> WEIGHT MODIFICATION
            else

              !first check to see if adjusted value not already calculated (skip if already calculated)
              if (valuep_padj(age,i,bignode,apnode,nodek)<-100.0) then
                !find weight for which valuep_married=valuep_divorced
                call find_valpwgt(valuep_divorced(age,ins_d,bignode,apnode,nodek),i,nodek)
              endif 
              valuep(age,i,bignode,m,apnode,nodek)=valuep_padj(age,i,bignode,apnode,nodek)
              valuek(age,i,bignode,m,apnode,nodek)=valuek_padj(age,i,bignode,apnode,nodek)
              consump(age,i,bignode,m,apnode,nodek)=conexpp_padj(age,i,bignode,apnode,nodek)
              consumk(age,i,bignode,m,apnode,nodek)=conexpk_padj(age,i,bignode,apnode,nodek)
              kappa(age,i,bignode,m,apnode,nodek)=kappa_padj(age,i,bignode,apnode,nodek)
              mu(age,i,bignode,m,apnode,nodek)=mu_padj(age,i,bignode,apnode,nodek)
              LFP(age,i,bignode,m,apnode,nodek)=LFP_padj(age,i,bignode,apnode,nodek)
              family(age,i,bignode,m,apnode,nodek)=famcare_padj(age,i,bignode,apnode,nodek)
              married(age,i,bignode,m,apnode,nodek)=12
              value(age,i,bignode,m,apnode,nodek)=mu(age,i,bignode,m,apnode,nodek)* &
                                                          valuep(age,i,bignode,m,apnode,nodek) + &
                                                        (1.0-mu(age,i,bignode,m,apnode,nodek))* &
                                                          valuek(age,i,bignode,m,apnode,nodek)

            end if
            end if !end if over limcom on or off
          enddo !end loop over insurance

          !insurance choices in period one:
		  if (ni==3) then !public option
			  if (age==1) then
				if (limcom==0) then
				  insurance(bignode,m,apnode,nodek)=insurance_married(bignode,m,apnode,nodek)
				else if (married(age,1,bignode,m,apnode,nodek)<15 .and. married(age,2,bignode,m,apnode,nodek)<15 .and. &
						 health.ne.1) then
				  insurance(bignode,m,apnode,nodek)=2
				else if (married(age,1,bignode,m,apnode,nodek)<15 .and. married(age,2,bignode,m,apnode,nodek)<15 .and. &
						 health==1) then
				  if (value(age,1,bignode,m,apnode,nodek)>value(age,2,bignode,m,apnode,nodek) .and. &
					  value(age,1,bignode,m,apnode,nodek)>value(age,3,bignode,m,apnode,nodek) ) then
					insurance(bignode,m,apnode,nodek)=1
				  elseif (value(age,3,bignode,m,apnode,nodek)>value(age,1,bignode,m,apnode,nodek) .and. &
					      value(age,3,bignode,m,apnode,nodek)>value(age,2,bignode,m,apnode,nodek) ) then
					insurance(bignode,m,apnode,nodek)=3 !public option
				  else 
					insurance(bignode,m,apnode,nodek)=2
				  endif
				else if (married(age,1,bignode,m,apnode,nodek)<15 .and. married(age,2,bignode,m,apnode,nodek)>15) then
				  insurance(bignode,m,apnode,nodek)=1
				else if (married(age,1,bignode,m,apnode,nodek)>15 .and. married(age,2,bignode,m,apnode,nodek)<15) then
				  insurance(bignode,m,apnode,nodek)=2
				else
				  insurance(bignode,m,apnode,nodek)=insurance_divorced(bignode,apnode,nodek)
				endif
			  endif
		  else	!not public option
			  if (age==1) then
				if (limcom==0) then
				  insurance(bignode,m,apnode,nodek)=insurance_married(bignode,m,apnode,nodek)
				else if (married(age,1,bignode,m,apnode,nodek)<15 .and. married(age,2,bignode,m,apnode,nodek)<15 .and. &
						 health.ne.1) then
				  insurance(bignode,m,apnode,nodek)=2
				else if (married(age,1,bignode,m,apnode,nodek)<15 .and. married(age,2,bignode,m,apnode,nodek)<15 .and. &
						 health==1) then
				  if (value(age,1,bignode,m,apnode,nodek)>value(age,2,bignode,m,apnode,nodek)) then
					insurance(bignode,m,apnode,nodek)=1
				  else
					insurance(bignode,m,apnode,nodek)=2
				  endif
				else if (married(age,1,bignode,m,apnode,nodek)<15 .and. married(age,2,bignode,m,apnode,nodek)>15) then
				  insurance(bignode,m,apnode,nodek)=1
				else if (married(age,1,bignode,m,apnode,nodek)>15 .and. married(age,2,bignode,m,apnode,nodek)<15) then
				  insurance(bignode,m,apnode,nodek)=2
				else
				  insurance(bignode,m,apnode,nodek)=insurance_divorced(bignode,apnode,nodek)
				endif
			  endif
		  endif
		  
        enddo !end loop over kid assets
      enddo !end loop over parent assets
    enddo !end loop over pareto weight         
    !**** END LOOP OVER PARETO WEIGHT *************************************************!

  enddo !end of loop over health (h)
  !*********************************************************************************************! 

  call MPI_ALLREDUCE(value(t,:,:,:,:,:),tempvec3,recsize3,MPI_DOUBLE_PRECISION,MPI_SUM,MPI_COMM_WORLD,ierr)
  value(t,:,:,:,:,:)=reshape(tempvec3,(/ni,nbignodes,nm,na,nonodesk/))
  call MPI_ALLREDUCE(valuep(t,:,:,:,:,:),tempvec3,recsize3,MPI_DOUBLE_PRECISION,MPI_SUM,MPI_COMM_WORLD,ierr)
  valuep(t,:,:,:,:,:)=reshape(tempvec3,(/ni,nbignodes,nm,na,nonodesk/))
  call MPI_ALLREDUCE(valuek(t,:,:,:,:,:),tempvec3,recsize3,MPI_DOUBLE_PRECISION,MPI_SUM,MPI_COMM_WORLD,ierr)
  valuek(t,:,:,:,:,:)=reshape(tempvec3,(/ni,nbignodes,nm,na,nonodesk/))

1234 continue !if skipped over cooperative solution

  call MPI_ALLREDUCE(valuep_divorced(t,:,:,:,:),tempvec4,recsize4,MPI_DOUBLE_PRECISION,MPI_SUM,MPI_COMM_WORLD,ierr)
  valuep_divorced(t,:,:,:,:)=reshape(tempvec4,(/ni,nbignodes,na,nonodesk/))
  call MPI_ALLREDUCE(valuek_divorced(t,:,:,:,:),tempvec4,recsize4,MPI_DOUBLE_PRECISION,MPI_SUM,MPI_COMM_WORLD,ierr)
  valuek_divorced(t,:,:,:,:)=reshape(tempvec4,(/ni,nbignodes,na,nonodesk/))

  call MPI_BARRIER(MPI_COMM_WORLD,ierr)
  
enddo !end of loop over time (t)
!***************************************************************************************************!

!gather all other matrices that will be necessary for simulations but were not necessary for solution

if (noncoop.ne.1) then

call MPI_ALLREDUCE(insurance,tempvec5,recsize5,MPI_DOUBLE_PRECISION,MPI_SUM,MPI_COMM_WORLD,ierr)
insurance=reshape(tempvec5,(/nbignodes,nm,na,nonodesk/))
call MPI_ALLREDUCE(insurance_divorced,tempvec6,recsize6,MPI_DOUBLE_PRECISION,MPI_SUM,MPI_COMM_WORLD,ierr)
insurance_divorced=reshape(tempvec6,(/nbignodes,na,nonodesk/))
call MPI_ALLREDUCE(mu,tempvec1,recsize1,MPI_DOUBLE_PRECISION,MPI_SUM,MPI_COMM_WORLD,ierr)
mu=reshape(tempvec1,(/nt,ni,nbignodes,nm,na,nonodesk/))

call MPI_ALLREDUCE(valuepNT0_married,tempvec1,recsize1,MPI_DOUBLE_PRECISION,MPI_SUM,MPI_COMM_WORLD,ierr)
valuepNT0_married=reshape(tempvec1,(/nt,ni,nbignodes,nm,na,nonodesk/))
call MPI_ALLREDUCE(valuepPT0_married,tempvec1,recsize1,MPI_DOUBLE_PRECISION,MPI_SUM,MPI_COMM_WORLD,ierr)
valuepPT0_married=reshape(tempvec1,(/nt,ni,nbignodes,nm,na,nonodesk/))
call MPI_ALLREDUCE(valuepFT0_married,tempvec1,recsize1,MPI_DOUBLE_PRECISION,MPI_SUM,MPI_COMM_WORLD,ierr)
valuepFT0_married=reshape(tempvec1,(/nt,ni,nbignodes,nm,na,nonodesk/))
call MPI_ALLREDUCE(valuepNT100_married,tempvec1,recsize1,MPI_DOUBLE_PRECISION,MPI_SUM,MPI_COMM_WORLD,ierr)
valuepNT100_married=reshape(tempvec1,(/nt,ni,nbignodes,nm,na,nonodesk/))
call MPI_ALLREDUCE(valuepPT100_married,tempvec1,recsize1,MPI_DOUBLE_PRECISION,MPI_SUM,MPI_COMM_WORLD,ierr)
valuepPT100_married=reshape(tempvec1,(/nt,ni,nbignodes,nm,na,nonodesk/))
call MPI_ALLREDUCE(valuepFT100_married,tempvec1,recsize1,MPI_DOUBLE_PRECISION,MPI_SUM,MPI_COMM_WORLD,ierr)
valuepFT100_married=reshape(tempvec1,(/nt,ni,nbignodes,nm,na,nonodesk/))

call MPI_ALLREDUCE(valuekNT0_married,tempvec1,recsize1,MPI_DOUBLE_PRECISION,MPI_SUM,MPI_COMM_WORLD,ierr)
valuekNT0_married=reshape(tempvec1,(/nt,ni,nbignodes,nm,na,nonodesk/))
call MPI_ALLREDUCE(valuekPT0_married,tempvec1,recsize1,MPI_DOUBLE_PRECISION,MPI_SUM,MPI_COMM_WORLD,ierr)
valuekPT0_married=reshape(tempvec1,(/nt,ni,nbignodes,nm,na,nonodesk/))
call MPI_ALLREDUCE(valuekFT0_married,tempvec1,recsize1,MPI_DOUBLE_PRECISION,MPI_SUM,MPI_COMM_WORLD,ierr)
valuekFT0_married=reshape(tempvec1,(/nt,ni,nbignodes,nm,na,nonodesk/))
call MPI_ALLREDUCE(valuekNT100_married,tempvec1,recsize1,MPI_DOUBLE_PRECISION,MPI_SUM,MPI_COMM_WORLD,ierr)
valuekNT100_married=reshape(tempvec1,(/nt,ni,nbignodes,nm,na,nonodesk/))
call MPI_ALLREDUCE(valuekPT100_married,tempvec1,recsize1,MPI_DOUBLE_PRECISION,MPI_SUM,MPI_COMM_WORLD,ierr)
valuekPT100_married=reshape(tempvec1,(/nt,ni,nbignodes,nm,na,nonodesk/))
call MPI_ALLREDUCE(valuekFT100_married,tempvec1,recsize1,MPI_DOUBLE_PRECISION,MPI_SUM,MPI_COMM_WORLD,ierr)
valuekFT100_married=reshape(tempvec1,(/nt,ni,nbignodes,nm,na,nonodesk/))

call MPI_ALLREDUCE(valueNT0_married,tempvec1,recsize1,MPI_DOUBLE_PRECISION,MPI_SUM,MPI_COMM_WORLD,ierr)
valueNT0_married=reshape(tempvec1,(/nt,ni,nbignodes,nm,na,nonodesk/))
call MPI_ALLREDUCE(valuePT0_married,tempvec1,recsize1,MPI_DOUBLE_PRECISION,MPI_SUM,MPI_COMM_WORLD,ierr)
valuePT0_married=reshape(tempvec1,(/nt,ni,nbignodes,nm,na,nonodesk/))
call MPI_ALLREDUCE(valueFT0_married,tempvec1,recsize1,MPI_DOUBLE_PRECISION,MPI_SUM,MPI_COMM_WORLD,ierr)
valueFT0_married=reshape(tempvec1,(/nt,ni,nbignodes,nm,na,nonodesk/))
call MPI_ALLREDUCE(valueNT100_married,tempvec1,recsize1,MPI_DOUBLE_PRECISION,MPI_SUM,MPI_COMM_WORLD,ierr)
valueNT100_married=reshape(tempvec1,(/nt,ni,nbignodes,nm,na,nonodesk/))
call MPI_ALLREDUCE(valuePT100_married,tempvec1,recsize1,MPI_DOUBLE_PRECISION,MPI_SUM,MPI_COMM_WORLD,ierr)
valuePT100_married=reshape(tempvec1,(/nt,ni,nbignodes,nm,na,nonodesk/))
call MPI_ALLREDUCE(valueFT100_married,tempvec1,recsize1,MPI_DOUBLE_PRECISION,MPI_SUM,MPI_COMM_WORLD,ierr)
valueFT100_married=reshape(tempvec1,(/nt,ni,nbignodes,nm,na,nonodesk/))

call MPI_ALLREDUCE(conexppNT0_married,tempvec1,recsize1,MPI_DOUBLE_PRECISION,MPI_SUM,MPI_COMM_WORLD,ierr)
conexppNT0_married=reshape(tempvec1,(/nt,ni,nbignodes,nm,na,nonodesk/))
call MPI_ALLREDUCE(conexppPT0_married,tempvec1,recsize1,MPI_DOUBLE_PRECISION,MPI_SUM,MPI_COMM_WORLD,ierr)
conexppPT0_married=reshape(tempvec1,(/nt,ni,nbignodes,nm,na,nonodesk/))
call MPI_ALLREDUCE(conexppFT0_married,tempvec1,recsize1,MPI_DOUBLE_PRECISION,MPI_SUM,MPI_COMM_WORLD,ierr)
conexppFT0_married=reshape(tempvec1,(/nt,ni,nbignodes,nm,na,nonodesk/))
call MPI_ALLREDUCE(conexppNT100_married,tempvec1,recsize1,MPI_DOUBLE_PRECISION,MPI_SUM,MPI_COMM_WORLD,ierr)
conexppNT100_married=reshape(tempvec1,(/nt,ni,nbignodes,nm,na,nonodesk/))
call MPI_ALLREDUCE(conexppPT100_married,tempvec1,recsize1,MPI_DOUBLE_PRECISION,MPI_SUM,MPI_COMM_WORLD,ierr)
conexppPT100_married=reshape(tempvec1,(/nt,ni,nbignodes,nm,na,nonodesk/))
call MPI_ALLREDUCE(conexppFT100_married,tempvec1,recsize1,MPI_DOUBLE_PRECISION,MPI_SUM,MPI_COMM_WORLD,ierr)
conexppFT100_married=reshape(tempvec1,(/nt,ni,nbignodes,nm,na,nonodesk/))

call MPI_ALLREDUCE(conexpkNT0_married,tempvec1,recsize1,MPI_DOUBLE_PRECISION,MPI_SUM,MPI_COMM_WORLD,ierr)
conexpkNT0_married=reshape(tempvec1,(/nt,ni,nbignodes,nm,na,nonodesk/))
call MPI_ALLREDUCE(conexpkPT0_married,tempvec1,recsize1,MPI_DOUBLE_PRECISION,MPI_SUM,MPI_COMM_WORLD,ierr)
conexpkPT0_married=reshape(tempvec1,(/nt,ni,nbignodes,nm,na,nonodesk/))
call MPI_ALLREDUCE(conexpkFT0_married,tempvec1,recsize1,MPI_DOUBLE_PRECISION,MPI_SUM,MPI_COMM_WORLD,ierr)
conexpkFT0_married=reshape(tempvec1,(/nt,ni,nbignodes,nm,na,nonodesk/))
call MPI_ALLREDUCE(conexpkNT100_married,tempvec1,recsize1,MPI_DOUBLE_PRECISION,MPI_SUM,MPI_COMM_WORLD,ierr)
conexpkNT100_married=reshape(tempvec1,(/nt,ni,nbignodes,nm,na,nonodesk/))
call MPI_ALLREDUCE(conexpkPT100_married,tempvec1,recsize1,MPI_DOUBLE_PRECISION,MPI_SUM,MPI_COMM_WORLD,ierr)
conexpkPT100_married=reshape(tempvec1,(/nt,ni,nbignodes,nm,na,nonodesk/))
call MPI_ALLREDUCE(conexpkFT100_married,tempvec1,recsize1,MPI_DOUBLE_PRECISION,MPI_SUM,MPI_COMM_WORLD,ierr)
conexpkFT100_married=reshape(tempvec1,(/nt,ni,nbignodes,nm,na,nonodesk/))

call MPI_ALLREDUCE(kappaNT0_married,tempvec1,recsize1,MPI_DOUBLE_PRECISION,MPI_SUM,MPI_COMM_WORLD,ierr)
kappaNT0_married=reshape(tempvec1,(/nt,ni,nbignodes,nm,na,nonodesk/))
call MPI_ALLREDUCE(kappaPT0_married,tempvec1,recsize1,MPI_DOUBLE_PRECISION,MPI_SUM,MPI_COMM_WORLD,ierr)
kappaPT0_married=reshape(tempvec1,(/nt,ni,nbignodes,nm,na,nonodesk/))
call MPI_ALLREDUCE(kappaFT0_married,tempvec1,recsize1,MPI_DOUBLE_PRECISION,MPI_SUM,MPI_COMM_WORLD,ierr)
kappaFT0_married=reshape(tempvec1,(/nt,ni,nbignodes,nm,na,nonodesk/))
call MPI_ALLREDUCE(kappaNT100_married,tempvec1,recsize1,MPI_DOUBLE_PRECISION,MPI_SUM,MPI_COMM_WORLD,ierr)
kappaNT100_married=reshape(tempvec1,(/nt,ni,nbignodes,nm,na,nonodesk/))
call MPI_ALLREDUCE(kappaPT100_married,tempvec1,recsize1,MPI_DOUBLE_PRECISION,MPI_SUM,MPI_COMM_WORLD,ierr)
kappaPT100_married=reshape(tempvec1,(/nt,ni,nbignodes,nm,na,nonodesk/))
call MPI_ALLREDUCE(kappaFT100_married,tempvec1,recsize1,MPI_DOUBLE_PRECISION,MPI_SUM,MPI_COMM_WORLD,ierr)
kappaFT100_married=reshape(tempvec1,(/nt,ni,nbignodes,nm,na,nonodesk/))

call MPI_ALLREDUCE(married,tempvec1,recsize1,MPI_DOUBLE_PRECISION,MPI_SUM,MPI_COMM_WORLD,ierr)
married=reshape(tempvec1,(/nt,ni,nbignodes,nm,na,nonodesk/))

endif

call MPI_ALLREDUCE(valuepNT_divorced,tempvec2,recsize2,MPI_DOUBLE_PRECISION,MPI_SUM,MPI_COMM_WORLD,ierr)
valuepNT_divorced=reshape(tempvec2,(/nt,ni,nbignodes,na,nonodesk/))
call MPI_ALLREDUCE(valuepPT_divorced,tempvec2,recsize2,MPI_DOUBLE_PRECISION,MPI_SUM,MPI_COMM_WORLD,ierr)
valuepPT_divorced=reshape(tempvec2,(/nt,ni,nbignodes,na,nonodesk/))
call MPI_ALLREDUCE(valuepFT_divorced,tempvec2,recsize2,MPI_DOUBLE_PRECISION,MPI_SUM,MPI_COMM_WORLD,ierr)
valuepFT_divorced=reshape(tempvec2,(/nt,ni,nbignodes,na,nonodesk/))
call MPI_ALLREDUCE(valuekNT_divorced,tempvec2,recsize2,MPI_DOUBLE_PRECISION,MPI_SUM,MPI_COMM_WORLD,ierr)
valuekNT_divorced=reshape(tempvec2,(/nt,ni,nbignodes,na,nonodesk/))
call MPI_ALLREDUCE(valuekPT_divorced,tempvec2,recsize2,MPI_DOUBLE_PRECISION,MPI_SUM,MPI_COMM_WORLD,ierr)
valuekPT_divorced=reshape(tempvec2,(/nt,ni,nbignodes,na,nonodesk/))
call MPI_ALLREDUCE(valuekFT_divorced,tempvec2,recsize2,MPI_DOUBLE_PRECISION,MPI_SUM,MPI_COMM_WORLD,ierr)
valuekFT_divorced=reshape(tempvec2,(/nt,ni,nbignodes,na,nonodesk/))

call MPI_ALLREDUCE(conexppNT_divorced,tempvec2,recsize2,MPI_DOUBLE_PRECISION,MPI_SUM,MPI_COMM_WORLD,ierr)
conexppNT_divorced=reshape(tempvec2,(/nt,ni,nbignodes,na,nonodesk/))
call MPI_ALLREDUCE(conexppPT_divorced,tempvec2,recsize2,MPI_DOUBLE_PRECISION,MPI_SUM,MPI_COMM_WORLD,ierr)
conexppPT_divorced=reshape(tempvec2,(/nt,ni,nbignodes,na,nonodesk/))
call MPI_ALLREDUCE(conexppFT_divorced,tempvec2,recsize2,MPI_DOUBLE_PRECISION,MPI_SUM,MPI_COMM_WORLD,ierr)
conexppFT_divorced=reshape(tempvec2,(/nt,ni,nbignodes,na,nonodesk/))
call MPI_ALLREDUCE(conexpkNT_divorced,tempvec2,recsize2,MPI_DOUBLE_PRECISION,MPI_SUM,MPI_COMM_WORLD,ierr)
conexpkNT_divorced=reshape(tempvec2,(/nt,ni,nbignodes,na,nonodesk/))
call MPI_ALLREDUCE(conexpkPT_divorced,tempvec2,recsize2,MPI_DOUBLE_PRECISION,MPI_SUM,MPI_COMM_WORLD,ierr)
conexpkPT_divorced=reshape(tempvec2,(/nt,ni,nbignodes,na,nonodesk/))
call MPI_ALLREDUCE(conexpkFT_divorced,tempvec2,recsize2,MPI_DOUBLE_PRECISION,MPI_SUM,MPI_COMM_WORLD,ierr)
conexpkFT_divorced=reshape(tempvec2,(/nt,ni,nbignodes,na,nonodesk/))


if (pid==0 .and. DoEstimation==0) then 

  open(unit=outcsv1,file='output/insurance.csv',action="write",status="replace")
  do m=1,nm
  do h=1,nbignodes
  do i=1,na
    write(outcsv1,"(100000(f0.12,',',:))") insurance(h,m,i,:)
  end do
  end do
  end do

  open(unit=outcsv1,file='output/insurance_divorced.csv',action="write",status="replace")
  do h=1,nbignodes
  do i=1,na
    write(outcsv1,"(100000(f0.12,',',:))") insurance_divorced(h,i,:)
  end do
  end do

  open(unit=outcsv1,file='output/valuep.csv',action="write",status="replace")
  do m=1,nm
  do h=1,nbignodes
  do i=1,na
    write(outcsv1,"(100000(f0.12,',',:))") valuep(1,1,h,m,i,:),valuep(1,2,h,m,i,:)
  end do
  end do
  end do
  
  open(unit=outcsv1,file='output/valuek.csv',action="write",status="replace")
  do m=1,nm
  do h=1,nbignodes
  do i=1,na
    write(outcsv1,"(100000(f0.12,',',:))") valuek(1,1,h,m,i,:),valuek(1,2,h,m,i,:)
  end do
  end do
  end do  

  open(unit=outcsv1,file='output/value.csv',action="write",status="replace")
  do m=1,nm
  do h=1,nbignodes
  do i=1,na
    write(outcsv1,"(100000(f0.12,',',:))") value(1,1,h,m,i,:),value(1,2,h,m,i,:)
  end do
  end do
  end do  

  open(unit=outcsv1,file='output/valuep_divorced.csv',action="write",status="replace")
  do h=1,nbignodes
  do i=1,na
    write(outcsv1,"(100000(f0.12,',',:))") valuep_divorced(15,1,h,i,:),valuep_divorced(15,2,h,i,:)
  end do
  end do

  open(unit=outcsv1,file='output/valuek_divorced.csv',action="write",status="replace")
  do h=1,nbignodes
  do i=1,na
    write(outcsv1,"(100000(f0.12,',',:))") valuek_divorced(15,1,h,i,:),valuek_divorced(15,2,h,i,:)
  end do
  end do
  
  open(unit=outcsv1,file='output/mu.csv',action="write",status="replace")
  do m=1,nm
  do h=1,nbignodes
  do i=1,na
    write(outcsv1,"(100000(f0.12,',',:))") mu(15,1,h,m,i,:),mu(15,2,h,m,i,:)
  end do
  end do
  end do  

  open(unit=outcsv1,file='output/married.csv',action="write",status="replace")
  do m=1,nm
  do h=1,nbignodes
  do i=1,na
    write(outcsv1,"(100000(f0.12,',',:))") married(15,1,h,m,i,:),married(15,2,h,m,i,:)
  end do
  end do
  end do  

  open(unit=outcsv1,file='output/assetp.csv',action="write",status="replace")
  do i=1,na
    write(outcsv1,"(100000(f0.12,',',:))") assetp(i,:)
  enddo
  open(unit=outcsv1,file='output/assetk.csv',action="write",status="replace")
  do i=1,na
    write(outcsv1,"(100000(f0.12,',',:))") assetk(i,:)
  enddo
endif

!*****************************************************!
!***********DONE SOLVING THE MODEL********************!
!*****************************************************!


!*****************************************************!
!***********SIMULATION********************************!
!*****************************************************!

print*, 'In simulation!'

valreturn=99999999
do m=1,nm

  call allocate_simulate(ier)
  call simmodel !simulation routine
  call writesim

  do momi=1,9
	moment_simu_vector(   momi)=wealthck50(1,momi)
	moment_simu_vector(9 +momi)=wealthck50(2,momi)
	moment_simu_vector(18+momi)=wealthck50(3,momi)	
  enddo
  do momi=1,7
  	moment_simu_vector(27+momi)=wealthck50(4,momi)
  enddo
  do momi=1,5
	moment_simu_vector(34+momi)=wealthck50(5,momi)
  enddo
  do momi=1,9
	moment_simu_vector(39+momi)=wealthcp50(1,momi)
	moment_simu_vector(48+momi)=wealthcp50(2,momi)
	moment_simu_vector(57+momi)=wealthcp50(3,momi)
  enddo 
  do momi=1,7
	moment_simu_vector(66+momi)=wealthcp50(4,momi)
  enddo
  do momi=1,5
	moment_simu_vector(73+momi)=wealthcp50(5,momi)	
  enddo
  moment_simu_vector(79)=avins(1)
  moment_simu_vector(80)=avins(2)
  moment_simu_vector(81)=avins(3)
  moment_simu_vector(82)=avins(4)
  moment_simu_vector(83)=avins(5)
  moment_simu_vector(84)=avinf(1)
  moment_simu_vector(85)=avinf(2)
  moment_simu_vector(86)=avinf(3)
  moment_simu_vector(87)=avinf(4)
  moment_simu_vector(88)=avinf(5)
  moment_simu_vector(89)=avSIpoverall
  moment_simu_vector(90)=avFTcohort(1)
  moment_simu_vector(91)=avFTcohort(2)
  moment_simu_vector(92)=avFTcohort(3)
  moment_simu_vector(93)=avFTcohort(4)
  moment_simu_vector(94)=avFTcohort(5)
    
  moment_err_vector=moment_data_vector-moment_simu_vector

  !Calculate objective function value (diagonally weighted)
  valreturnmat(m)=0.d0
  do momi=1,Nmoments
    valreturnmat(m)=valreturnmat(m)+weight_matrix(momi,momi)*moment_err_vector(momi)**2
  enddo

  if (pid==0) then
    print*,'value of this m',musimstartparam,valreturnmat(m)
    print*,'ins dist',avins(:)
    print*,'ins PO dist',avinsPO(:)
    print*,'inf dist',avinf(:)
    print*,'medicaid',avSIpoverall
	print*,'full time / part time', avFToverall, avPToverall
    print*,'parent assets avg',wealthcp50(1,:)
    print*,'parent assets avg',wealthcp50(5,:)
    print*,'child assets avg',wealthck50(1,:)
    print*,'child assets avg',wealthck50(5,:)
  endif

  goto 1111 !skip other values of m
enddo !loop over musimstartparam values


  !*****************************************************!
  !***********END SIMULATION****************************!
  !*****************************************************!

1111 continue

  call system_clock ( t2, clock_rate, clock_max )
  if (pid==0) write(*,*) 'Done! Elapsed real time = ', real ( t2 - t1 ) / real ( clock_rate )

  estmatrix(whichiter,1)=formalpref
  estmatrix(whichiter,2)=stigma
  estmatrix(whichiter,3)=altruism
  estmatrix(whichiter,4)=alpha
  estmatrix(whichiter,5)=Lshifter
  estmatrix(whichiter,6)=gammab
  estmatrix(whichiter,7)=gamma
  estmatrix(whichiter,8)=cminp
  estmatrix(whichiter,9)=mcaidpref
  estmatrix(whichiter,10)=valreturn

  if (pid==0) then
    do i=1,whichiter
      print '(f18.8 f18.8 f18.8 f18.8 f18.8 f18.8 f18.8 f18.8 f18.8)', estmatrix(i,:)
    enddo
  endif
  whichiter=whichiter+1

end subroutine runmyfunction